Current-Phase Relation of a Bose-Einstein Condensate 
Flowing Through a Weak Link 
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We study the current-phase relation of a Bose-Einstein condensate flowing through a repulsive 
square barrier by solving analytically the one dimensional Gross-Pitaevskii equation. The barrier 
height and width fix the current-phase relation j(5(f>), which tends to j ~ cos(<50/2) for weak barriers 
and to the Josephson sinusoidal relation j ~ sin(<5</>) for strong barriers. Between these two limits, 
the current-phase relation depends on the barrier width. In particular, for wide enough barriers, 
we observe two families of multivalued current-phase relations. Diagrams belonging to the first 
family, already known in the literature, can have two different positive values of the current at 
the same phase difference. The second family, new to our knowledge, can instead allow for three 
different positive currents still corresponding to the same phase difference. Finally, we show that the 
multivalued behavior arises from the competition between hydrodynamic and nonlinear-dispersive 
components of the flow, the latter due to the presence of a soliton inside the barrier region. 

PACS numbers: 03.75.Lm, 74.50.+r, 47.37.+q 



I. INTRODUCTION 

The current-phase relation characterizes the flow of a su- 
perfluid/superconductor through a weak link [H, @, 
The latter is a constricted flow configuration that can be 
realized in different ways: i) apertures in impenetrable 
walls mostly for helium, ii) sandwich or bridge structures 
for superconductors, and iii) penetrable barriers created 
by laser beams for ultracold dilute gases. Much infor- 
mation about such systems can be extracted from the 
current-phase relation, which, given a fluid, depends only 
on the link properties. For instance, with 4 He, the tran- 
sition from the usual AC Josephson effect to a quantized 
phase slippage regime 0] corresponds to the switching 
from a sine-like current phase relation to a multivalued 
one [f|. 

A weak link configuration can be modelled very generally 
upon taking a portion of a superfluid/supcrconductor to 
have "different conduction properties" with respect to 
the rest of the system. Two pieces of superconductor 
joined by a third superconducting region with a smaller 
coherence length provide one example, whose current- 
phase relation in one dimension has been studied with 
the Ginzburg-Landau equation Q. 

In the context of ultracold dilute gases, raising a repul- 
sive penetrable barrier across the flow yields an equiva- 
lent configuration. For instance, with Bose-Einstein con- 
densates (BEC), Josephson effect(s) have been theoreti- 
cally studied 0,1, i, liol , lTi\ and experimentally demon- 
strated using multiple well traps [l2|, 0, EH ■ Theoret- 
ically, the current-phase relation has been studied for a 
flow through a repulsive square well with fermions across 
the BCS-BEC crossover by means of one dimensional 
Bogoliubov-de Gennes equations [IH, for weak barriers 
with bosons in a local density approximation (l6| . and 
for fermions on the BEC side of the crossover using a 
nonlinear Schrodinger equation approach [13]. 
In this manuscript, wc study the current-phase relation 
for a BEC flowing through a repulsive square well. The 



weak link configuration, and in turn the current-phase 
relation, is then determined by the barrier height with 
respect to the chemical potential and by the barrier 
width with respect to the healing length. Though we 
solve a one-dimensional Gross-Pitaevskii equation, the 
results presented in this manuscript are not just rele- 
vant for BECs, but also include the essential features of 
current-phase relations of superconducting or superfluid 
He-based weak links when governed by the Ginzburg- 
Landau equation. For any barrier width, we find that in 
the limit of zero barrier height, the current phase relation 
tends to j(S(f>) — cos(<50/2), with c,^ being the bulk 
sound velocity, which corresponds to the phase across a 
grey soliton at rest with respect to the barrier. On the 
other hand, if the barrier height is above the bulk chem- 
ical potential at zero current, the limit of tunneling flow 
is reached either when the barrier height is much bigger 
than the bulk chemical potential at zero current or when 
the barrier width is much larger than the bulk healing 
length. In this regime, we recover the the usual Joseph- 
son sinusoidal current-phase relation and obtain an ana- 
lytical expression for the Josephson critical current as a 
function of the weak link parameters. For barriers wider 
than the healing lenght inside the barrier region, we ob- 
serve two families of multivalued (often called reentrant) 
current-phase relations. The first, already studied since 
the early works on superconductivity 0, [y] , shows a posi- 
tive slope of the current when the phase difference is close 
to 7T, thereby reaching a phase difference larger than tt 
at least for small currents. The second family, appear- 
ing at a smaller barrier height, has instead a negative 
slope of the current close to 7T, and in some cases can 
remain within the — tt interval across the whole range 
of currents. These two families can also be distinguished 
by the maximum number of different positive currents 
corresponding to the same phase difference: two for the 
first family, three for the second one. As the first kind 
of reentrant behavior was proven to be connected to the 
onset of phase-slippage in the AC Josephson effect 
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the second might then be connected to the appearance 
of new features in the Josephson dynamics. We finally 
observe that the hysteresis characterizing both families 
of reentrant current-phase relations is always due to the 
competition between a hydrodynamic component of the 
flow and a nonlinear-dispersive component, the latter due 
to the presence of a soliton inside the barrier region. The 
two components can coexist only for barriers wide enough 
to accomodate a soliton inside. In this spirit, we develop 
a simple analytical model which describes very well reen- 
trant regimes of current-phase relations. 



II. THEORETICAL DEVELOPMENTS 

A. Stationary Solutions 

We consider a dilute repulsive Bose-Einstein condensate 
at zero temperature flowing through a ID rectangular 
potential barrier. We look for stationary solutions of the 
ID GPE 0: 



2m 



d xx y + V ext (x)^ + g\y\ 2 y = M *, (1) 



where vE'(x) = \f n(x) exp[i(j)(x)] is the complex order 
parameter of the condensate, [i is the chemical poten- 
tial, and g = 4irh 2 a s /m with m the atom mass and 
a s > the s-wavc scattering length. The order param- 
eter phase 4>{x) is related to the superfluid velocity via 
v(x) = (h/m)d x (f>(x). The piecewise constant external 
potential describes the rectangular barrier of width 2d 
and height Vq: 



V cxt (x) = V , \x\ < d , 
V cxt (x) = , \x\ > d . 



(2) 



We consider solutions of Eq. (fTJ) which are symmetric 
with respect to the point x = 0, therefore discarding cases 
in which a reflected wave is present . Such symmetric 
solutions in the presence of a barrier exist due to the 
nonlincarity in the GPE. We also restrict our analysis to 
subsonic flows Voo < Coo , with = ^Jgn^/m being the 
sound velocity for a uniform condensate of density n^. 
As boundary conditions, we fix the condensate density 
rioo and velocity at x = ±oo, thereby determining the 
chemical potential fi = gn^ + ^mv 2 ^. Using the relation 
^ = Y / ncxp[i0], Eq. |T]) can be split into a continuity 
equation, stating the constancy in space of the current 
j = n(x)v(x) = TiaoVco, and an equation for the density 
only 
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where we have used the continuity equation v(x) = 
j /n(x). Its solution n(x) is expressed in terms of Jacobi 
elliptic functions 0, l2l| . Due to symmetry, we need only 
consider half of the space x > 0. The solution outside 
the barrier x > d becomes 



n out {x) = ?ioo-Ax>+A>otanh 2 



mg 

-^-Aoo(x - d) + xo 



where xq = arctanh-^/ (rid — / g) / '^oo , ^4oo = ^oo — 
mv to/9 ^ Oj and rid is the density at the barrier edge 
x = d. The solution inside the barrier x < d is: 



rii(x) = no + A 



n 2 {x) = n + A; 



sn 2 [bix, ki] 
cn 2 [6ia;, fci] 
1 - cn[6 2 a;, k 2 ] 
1 + cn[b 2 x, k 2 ] 



A > and A x > , 



else 



(5) 

where uq is the density at x = 0, A = (no — 2p,/g) 2 — 



Amj / grig with fx = fj, — Vo, A± = 3no/2 — jl/g — vA/2, 
and A 2 = ^/2(ng — n jl/g + mj 2 /2gn ). Finally, the pa- 
rameters entering the Jacobi sines sn and cosines cn are 

&i = y/mgiVA + A^/h 2 , h = (VA/ '(VA+AO) 1 / 2 > and 
b 2 = ^jAmgA 2 /h\ k 2 = {{A 2 -A,- A /A/2)/2A 2 ) 1 / 2 . 
Given and v^, we are left with two free parame- 
ters: rig and rid, which are next determined by matching 
the density and its derivative at the barrier edge x = d. 
First, the derivative matching condition, enforced using 
the first integral of Eq. ([3]), allows us to write rid as a 
function of no for any value of A: 
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(4) 



(6) 

The density matching equation rid — ni(d) (i = 1,2) is 
then solved by a numerical root finding method, yelding 
no. 

Two bounded solutions are always found [22| , an example 
is given in the upper panel of Fig. [T] In the following, 
the solution which tends to a plane wave for Vo — > 
will be referred to as "upper solution" , while the one 
tending to a grey soliton will be called "lower solution" 
[20j | . For given barrier parameters Vq and d, and at a 
fixed density at infinity noo, the solutions exist up to a 
critical injected velocity Voq = v c < Coo, at which they 
merge and disappear. This behavior was found in the 
case of a repulsive delta potential in [23[. Similarly, in 
fl9l ]. the same kind of merging was reported for a ID BEC 
flow through a repulsive square well when the width of 
the latter increases. 



B. Current-phase Relation 

As pointed out in the Introduction, the current-phase 
relation j(S(f>) for a given superfluid only depends on the 
properties of the weak link, in our case the barrier height 
Vo with respect to the chemical potential and the width 
2d with respect to the healing length. For a fixed current 
j = v(x)n(x), the phase difference across the system is 
calculated using the relation <j)(x) = J dy(m/h)j /n(y) 
and then rcnormalizcd by the phase accumulated by the a 
plane wave with the same boundary conditions in absence 
of barrier (see lower panel in Fig. [1]) . Two different values 
of 5<p are found, corresponding to the upper and lower 
solutions. 

In this section, we will use dimensionless quantities, 
employing the chemical potential at zero current grioo 
as the unit of energy, the bulk healing length ^ = 
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VoMoo = 0.4 
3 1 'curiae = 0.3 

— Upper solution 

— Lower solution 



5<f> 



FIG. 1: (Color online) Typical solutions for a barrier 
with width 2d = 4^. Density (upper panel) and phase 
(lower panel) as a function of position are shown for both 
the upper and the lower solution (see inset). Dashed lines 
in the lower panel correspond to the phase accumulated 
by a plane wave in absence of barrier. 



h/y / 2mgn 00 as the unit of length, and h/g tiqq as the 
unit of time. Exploiting the symmetry of the system 
about x = 0, the phase difference can be written as Sfa = 
lim x ^ oo [J jdy/n i (y) + J^jdy/n out (y)-jx/n 00 ] ,i = 1,2, 
where the third term is the phase difference accumulated 
by the plane wave. The limit can be calculated using 
Eq. (@|, yielding 



jdy 

m(y) 



-jd+2 



arcsin( 



J 



J 



(7) 

The first two terms in Eq. ([7]) correspond to the phase 
acquired inside the barrier while the third, which we can 
call the pre-bulk term, gives the phase accumulated out- 
side the barrier, where the density has not yet reached 
its bulk value rioo. We have taken the latter to be one. 
In Fig. [21 the current phase relation is shown for different 
barrier widths and heights. Each curve has a maximum 
at the point (Sfa,j c ), with j c = n^Vc being the critical 
current at which the two stationary solutions merge and 
disappear. The upper solutions constitute the part of 
the current-phase diagram which connects the maximum 
with the point (8<j> = 0,j — 0), while the lower ones 
belong to the branch connecting the maximum to the 
point (6<j) = 7r,j = 0). Indeed, we will now show that, 
for any d and Vo — ► 0, the upper branch tends to a plane 
wave, while the lower branch tends to a grey soliton. In 
order to have a finite n^ in this limit, the term in square 
brackets in Eq. (JS|) must tend to zero, yielding a cubic 
equation with two coincident solutions no = 1 and a third 
no = j 2 /2, where we have set noo = 1 for simplicity. For 
?io = 1 we have A = (1 — j 2 /2) 2 > 0, corresponding 
to the plane wave solution n\{x) = uq = 1. For uq = 
j 2 /2 we obtain instead A = 0, A\ = j 2 /2 - 1 < 0, 
corresponding to a grey soliton solution ri2(x) = n$ + 
(1 — no) tanh 2 (y / l/2 — no/2x). Therefore, in this limit 



d/£oc = 


= 0.5 


■ Vo/gn^ 


= 0.001 


♦ V /gn„ 


= 0.2 


• Vo/gn x 


= 0.7 


A Vo/S"o= 


= 16 



3_ 





= 3 


■ Vo/s»=c 


= 0.001 


♦ V /gn x 


= 0.035 


• VJ/jTfe, 


= 0.065 




= 1.8 




jr/2 3?r/4 



FIG. 2: (Color online) Current-phase relation for a bar- 
rier with half width d = 0.5£oo (left panel) and d = 3£oo 
(right panel), for different barrier heights Vq (sec inset). 
The arrows sketch the behavior of the maximum of the 
curves upon increasing Vq. 



5<f>i = for any j, meaning that the upper branch is 
actually a vertical line, while for the lower branch we 
have 



cos( — ) = — — . 



(8) 



This curve has a maximum at 5 fa = 0, corresponding to 
j = s/2, that is, the sound velocity Coo in dimensionless 
units. 

The arrows in Fig.[5]skctch the behavior of the maximum 
of the current-phase relation (8<p Cl j c ) as the height Vo is 
increased at a fixed barrier width 2d. For any width, the 
current-phase diagram initially takes a cosine-like shape 
(Eq. {8}) when Vb ~ (red squares in Fig. [2|), and tends 
to a sin((5</>) shape for sufficiently large Vo, characterizing 
the Josephson regime of tunneling flow (blue triangles 
in Fig. [2|). Between these two limits, the behavior of 
the maximum is determined by the barrier width. For 
thin barriers (d < £oo), as shown in the left panel of 
Fig. [21 the point (5fa,j c ) moves down-right, reaching 
the Josephson regime keeping 8<j) c always smaller than 
7r/2. On the other hand, for sufficiently wider barriers, 
5 fa is able to reach values larger than tt/2 within a fi- 
nite range of Vo during the down-right displacement of 
the maximum. The latter then moves down-left to finally 
enter the Josephson regime, as shown in the right panel 
of Fig. [2l We note that in this way, while Vq is increased, 
the phase 8 fa takes the value 717 '2 twice, but only the sec- 
ond time entering the Josephson regime with a sinusoidal 
current-phase relation. The first time (orange diamonds 
in the right panel of Fig. [2]) the flow is not yet in the tun- 
neling regime since Vo is much smaller than the chemical 
potential, indeed the current-phase relation is symmetri- 
cal with respect to tt/2, but not sinusoidal. 
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C. Josephson Regime 

As described in the previous section, for strong enough 
barriers the flow enters the tunneling regime, and the 
current-phase relation takes a sinusoidal form. In the 
following, we will describe analytically this behavior, de- 
riving a relation between the Josephson critical current 
and the barrier parameters Vq and d. Since we are now 
interested in tunneling flows, we will take Vq > gn^ and 
will show that the Josephson regime is attained by either 
increasing the barrier height Vq or its width 2d. 
Since in this regime the injected velocity of a stationary 
flow must be much smaller than the sound veloc- 
ity Coo, the chemical potential can be written as fi ~ 
grioo — Vq < 0, upon neglecting the kinetic energy term 
mv 2 <J /2. Moreover, the density inside the barrier no being 
exponentially small, we can write A ~ (2(l/g) 2 — s > 0, 
with s = 4(no// + mj 2 /no)/g, where we have neglected 
Hq. Thus, both the upper and lower solutions are of the 
kind 774(2;), with A\ ~ uq — mj 2 /2riQfl, b\ ~ ^J2m\jl\/h 2 , 
and k ~ 1. The density in the Josephson regime has thus 
the form 



77 jos (x) = n + (n + - — — 7) sum 
2n /-t 




(9) 



In order to write the density matching equation = 
774 (d), we approximate by discarding both n\ and 
2mv 2 n 00 /g in Eq. ([S]), obtaining a quadratic equation 

for 77,0. Further assuming that sinh 2 -i/— 3> 1, with 

€d = h 2 /2md being the kinetic energy associated with 
the barrier length scale d, the solutions of the above equa- 
tion are of the form n r 
2 



4'- 



7io(l ± \/l — q) with 7io = 
(5^00/4^0)7700/ sinh 2 y/\p,\/ed, and q = mj 2 /2nl\fl\. 
The Josephson critical current corresponds to the merg- 
ing of the two solutions at q = 1, and reads 



Jjos 



2|^| gn<x 
m Vq 



-2y/\ji\/e 



(10) 



where we have used sinh(x) ~ exp(x)/2, for x ^> 1. The 
critical velocity for a bosonic and fermionic supcrfiuid 
flowing through a repulsive square well has been calcu- 
lated in [l6j within the local density approximation (for 
BEC case see also [HI)- Analytical expressions for the 
critical current of a BEC flow were found for both slowly 
varying and weak barriers [23[. Eq. pTT)) thus enriches 
the above set of analytical results by providing the criti- 
cal current for strong barriers. 

Finally, we calculate the current-phase relation using 
Eq. ([7])- Since jj OS is exponentially small, only the first 
term in Eq. ([7]) contributes, and the integral can be per- 
formed analytically to obtain: 



J_ 



7r — arcsm 



_3_ 

J J'J, 



(11) 



Thus, we recover the sinusoidal current-phase relation 
characterizing a Josephson regime of tunneling flow. 
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FIG. 3: (Color online) Current-phase relation in the 
reentrant regime d = 20£oo- Insets show in more detail 
the shape of the diagram close to the critical point. The 
arrows sketch the behavior of the maximum of the curves 
upon increasing Vq. 
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FIG. 4: (Color online) Typical density profiles in the 
reentrant regime d = 20£oo for the upper and lower solu- 
tions. 

D. Reentrant Regimes 

When the barrier width 2d greatly exceeds the healing 
length, the current-phase relation becomes multivalued 
as shown in Fig. [3] for d = 20£oo ■ These so-called reen- 
trant current-phase diagrams were first predicted for long 
superconducting weak links @,@|. Remarkably, they have 
been experimentally demonstrated with superfluid 4 He 

In our system, two kinds of multivalued diagrams are 
found for a fixed barrier width 2d. We designate the first 
kind (green dots), appearing at larger barrier heights Vq, 
as reentrant type I, and the second kind (orange dia- 
monds in Fig. [3]), occurring for smaller Vo, as reentrant 
type II. They differ in the behavior of the lower branch 
at small currents j <C 1. The phase decreases with in- 
creasing j in type II diagrams while it increases in type 
I, reaching values larger than n for j ^ 1. These two 
families of diagrams also differ in the number of positive 
values of the current j corresponding to the same phase 
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difference S(f>. Indeed, diagrams of type I can have two 
values of j at the same 5(j), while diagrams of type II 
can allow for three. The existence of these type II dia- 
grams, to our knowledge, has not been discussed in the 
literature. In only type I current-phase relations were 
observed and connected to the onset of phase-slip dissi- 
pation in the system. In this spirit, a type II diagram 
might for instance lead to an instability of a different 
kind. 

In the remainder of this section, we will develop an an- 
alytical model that captures the essential features un- 
derlying both kinds of reentrant current-phase diagrams. 
Examination of typical density profiles belonging to the 
reentrant regime (see Fig. @| suggests to construct the 
lower solution by starting from the upper one at the 
same current, then adding a grey soliton inside the bar- 
rier region. Since we are dealing with wide barriers, we 
describe the upper solution in the local density approxi- 
mation (LDA) (its current-phase relation in this approx- 
imation is also discussed in [HI). At a fixed current j, 
the density inside the barrier \x\ < d is constant and 
equal to no = j±/Zg + 2p,cos(ui/3)/3g, when j < j t h, or 
no = fl/3g + 2/icos(7r/3 — w/3)/3<7, when j > j t h with 
jth defined by j 2 h = 4fl 3 /27mg 2 , and to = arccos |1 — 
27mg 2 j 2 /4/i 3 |. The phase difference calculated within 
LDA misses the pre-bulk term in Eq. ((7]), thus, for the 
upper branch, it is simply 5<f>i = 2mjd(l/no — l/n^/h. 
Now, for the density profile of the lower branch 712(21), 
we take a grey soliton (Eq. Q) placed inside the barrier 
at x = 0, with a bulk density given by no and a bulk 
velocity vq = j /no while in the region |x| > d we keep 
the density profile of the upper branch, that is, a con- 
stant density rioo and velocity j/n^. Notice that in this 
section no stands for the density of the upper solution at 
x = 0, as indicated in Fig. [4] The density at x = for 
the lower solution corresponds to the center of the dip in 
the grey soliton density profile. Finally, using Eq. (J7]) we 
obtain the phase difference for the lower branch 



2 arccos i 



gn 2 (d) 



(12) 



At a given current j, the overall phase difference cor- 
responding to the lower solution has two contribu- 
tions: 1) the "hydrodynamics phase" 8<pi coming from 
LDA and 2) the "nonlinear-dispersive phase" S<fi so \ = 
2 arccos y/ mv 2 / gn 2 (d) accumulated across the grey soli- 
ton. While 6(f>i is a monotonically increasing function 
of j, 5<f) so \ is instead monotonically decreasing, starting 
from 7r at zero current [25| . 

Therefore, the hysteresis characterizing a reentrant 
current-phase relation is due to the competition between 
the hydrodynamic and the nonlinear-dispersive compo- 
nents of the flow, which can coexist only for barriers 
wide enough to accomodate a soliton inside. In par- 
ticular, we can derive a condition for the appearance 
of type I reentrant behavior upon expanding Eq. (I12|) 
for small currents, and requiring 5<j) 2 > tt- Using 
arccos(a;) ~ ir/2 — x, for a; < 1, we get 8<p 2 — tt + 2jis, 
where v = md(l/rio — l/ n °o)/^ — V m l 9 n %-> an d wc 
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FIG. 5: (Color online) Comparison between the 
LDA+soliton model (left panel) or the Deaver-Pierce 
model (right panel) and the exact results. 



have taken n 2 (d) ~ no. The condition for type I reen- 
trance to appear is thus v > 0. For j « 1, we have 
no — [ijg — n x — Vo/g, and since within the LDA ap- 
proximation Vo *C grim, the condition v > takes the 
simple form 



Vo d 



> 1, 



(13) 



with £o = hj \J2mgno being the healing length inside the 
barrier region where the density is no. Equation Q13p . 
holding for Vq <C gn^,, has a clear physical meaning: 
in order to have a type I reentrant current phase dia- 
gram, the barrier width 2d must be sufficiently larger 
than 2"\/2£o, which is the characteristic size of a soliton 
placed inside the barrier. 

In the left panel of Fig. [SJ we compare the current-phase 
relation calculated with the above model (solid lines) to 
the exact results. Within the reentrant regime, for both 
type I and type II, the agreement is striking with only 
slight differences close to the the maximum (6(j) c ,j c ). On 
the other hand, for thin/strong barriers, LDA, and in 
turn the above model, is in clear disagreement with the 
exact results. (Sec cases d = 3£oo in the left panel of 
Fig.©. 

The difference between type I and type II diagrams has a 
physical interpretation within the above model, namely 
that the hydrodynamic component of the flow dominates 
for all currents in type I reentrance (excluding the re- 
gion j ~ j c [26|), while it is overcome by the nonlinear- 
dispersive part for sufficiently small current in type II. 

In the literature 0, multivalued current-phase relations 
are typically modelled by describing the weak-link with 
an equivalent circuit containing a linear inductance in se- 
ries with a sinusoidal inductance, the latter correspond- 
ing to an ideal Joscphson junction. When compared to 
our GPE exact results (see right panel in Fig. [5]), this 
model in general fails to describe the curvature of both 
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branches, and misses type II rcentrance, as well as nearly 
free regimes (e.g. red squares in Fig. I2I3|) since it does not 
allow the phase 54> c , corresponding to the maximum of 
the diagram, to be smaller than tt/2. It proves quite ac- 
curate only for sufficiently large barrier heights Vq, very 
close to the Josephson regime. In the 4 He experiment 
of 0] , this so called Dcaver-Pierce model [27[ agrees well 
with the measured current-phase relation. This might 
be due to the fact that these experiments are performed 
with a fixed weak link configuration, moving the system 
across the transition between Josephson and type I reen- 
trant regime, upon changing the He healing length with 
temperature, but always staying sufficiently close to the 
Josephson regime. 

III. CONCLUSIONS 

We studied the current-phase relation of a BEC flow- 
ing through a weak link created by a repulsive barrier. 
The link is thus modelled by two parameters, the barrier 
height and width, which fix the current-phase relation. 
Even though we solved a simplified model, wc believe 
that the results obtained will also be relevant for super- 
conducting and superfluid He-based weak links. 
We obtained analytical results for the weak barrier limit, 
for which the current-phase relation has a ~ cos(<50/2) 
form, and for the strong barrier limit, for which it takes 



a ~ sin(<5</>) form characterizing the Josephson regime. 
In particular, we derived an expression for the Josephson 
critical current as a function of the link parameters. 
Finally, we found two kinds of multivalued current-phase 
diagrams which we show, by means of an analytical 
model, to appear due to the competition between a hy- 
drodynamic and a nonlinear-dispersive part of the flow, 
which can coexist only for barriers wide enough to ac- 
comodate a soliton inside. The first kind, showing two 
different positive currents at same phase difference, has 
recently been experimentally demonstrated with 4 He and 
proven to be connected to the appearance of phase slip- 
page in the AC Josephson dynamics. The second one, 
new to our knowledge, can instead allow for three dif- 
ferent positive currents corresponding to the same phase 
difference. We believe that this new kind of hysteresis 
in the current-phase relation can be associated with new 
features emerging in the Josephson dynamics, which will 
be studied elsewhere. 
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